
restore,'./masses/mfn_mega_td1.sav',/ver

help,dist,/str

nreal = n_elements(dist.mass[0])
print,nreal

cgPlot,[0],[0],/nodata,xr=[0,5],yr=[1d-1,50],xtit='Log!d10!n M',ytit='N',$
       /ylog,/yst,ytickformat='exponent10'


hist = fltarr(nreal,51)
;; For each realization, make distance cuts, make histogram
FOR i=0,nreal-1 DO BEGIN
   
   d = dist[i].objs / 1.d3
   dind = where(d GE 1. AND d LE 8., nd)
   
   hist[i,*] = histogram(alog10(dist[i].mass[dind]),binsize=0.1d,$
                                min=0.00,max=5.0,loc=xh)
   xp = xh + 0.5*(median(xh[1:*]-xh[0:*]))                         
   ;; cgOplot,xp,hist[i,*],psym=10
ENDFOR


;; Loop through bins to get mu, sig
muhist = fltarr(51)
sighist = fltarr(51)

FOR j=0,51-1 DO BEGIN
   
   muhist[j] = mean(hist[*,j])
   sighist[j] = stddev(hist[*,j])
   
ENDFOR
yf = mpfitpeak(xp,muhist,A,nterms=3)
cgOplot,xp,yf,color='red',thick=3


vline,alog10(80),thick=3,color='lime green',/log
vline,alog10(7d3),thick=3,color='lime green',/log
vline,/h,0.5,thick=3,linestyle=3,color='deep pink'
;; vline,/h,1.0,thick=3,linestyle=3,color='deep pink'

cgOplot,xp,muhist,psym=10
oploterror,xp,muhist,sighist,errcolor='dodger blue',psym=10


;; plothist,alog10(mfn),xarr,yarr,bin=0.1

yh = histogram(alog10(mfn),binsize=0.1d,min=0.00,max=5.0,loc=xh)
xp = xh + 0.5*(median(xh[1:*]-xh[0:*]))                         
;; cgOplot,xp,yh,psym=-16,color='deep pink'        

print,n_elements(xp)


END
